

###################################################################

  Table 5. Comparisons of Test Statistics and Fit Indices

###################################################################


library(lavaan)

ess<-read.csv("/Users/bv533/Dropbox/LM-Test/Empirical Application/ESS/ESS.csv")

ess<-as.data.frame(ess)

DE=ess[ess$cntry=="DE",]


####################

  Cleaning data

###################


#1 Equality

DE$ipeqopt
DE$ipeqopt[DE$ipeqopt==7]<-NA
DE$ipeqopt[DE$ipeqopt==8]<-NA
DE$ipeqopt[DE$ipeqopt==9]<-NA

#2 Understand

DE$ipudrst
DE$ipudrst[DE$ipudrst==7]<-NA
DE$ipudrst[DE$ipudrst==8]<-NA
DE$ipudrst[DE$ipudrst==9]<-NA


#3 Environment

DE$impenv
DE$impenv[DE$impenv==7]<-NA
DE$impenv[DE$impenv==8]<-NA
DE$impenv[DE$impenv==9]<-NA

#4 Help others

DE$iphlppl
DE$iphlppl[DE$iphlppl==7]<-NA
DE$iphlppl[DE$iphlppl==8]<-NA
DE$iphlppl[DE$iphlppl==9]<-NA


#5 Loyalty

DE$iplylfr
DE$iplylfr[DE$iplylfr==7]<-NA
DE$iplylfr[DE$iplylfr==8]<-NA
DE$iplylfr[DE$iplylfr==9]<-NA



#6 Modesty
DE$ipmodst
DE$ipmodst[DE$ipmodst==7]<-NA
DE$ipmodst[DE$ipmodst==8]<-NA
DE$ipmodst[DE$ipmodst==9]<-NA


#7 Tradition
DE$imptrad
DE$imptrad[DE$imptrad==7]<-NA
DE$imptrad[DE$imptrad==8]<-NA
DE$imptrad[DE$imptrad==9]<-NA


#8 Rules
DE$ipfrule
DE$ipfrule[DE$ipfrule==7]<-NA
DE$ipfrule[DE$ipfrule==8]<-NA
DE$ipfrule[DE$ipfrule==9]<-NA


#9 Propriety

DE$ipbhprp
DE$ipbhprp[DE$ipbhprp==7]<-NA
DE$ipbhprp[DE$ipbhprp==8]<-NA
DE$ipbhprp[DE$ipbhprp==9]<-NA


#10 Security
DE$impsafe
DE$impsafe[DE$impsafe == 7]<-NA
DE$impsafe[DE$impsafe == 8]<-NA
DE$impsafe[DE$impsafe == 9]<-NA

#11 Strong Government

DE$ipstrgv
DE$ipstrgv[DE$ipstrgv==7]<-NA
DE$ipstrgv[DE$ipstrgv==8]<-NA
DE$ipstrgv[DE$ipstrgv==9]<-NA


#12 Race
DE$imdfetn
DE$imdfetn[DE$imdfetn ==7]<-NA
DE$imdfetn[DE$imdfetn ==8]<-NA
DE$imdfetn[DE$imdfetn ==9]<-NA


#13 EU poor

DE$eimpcnt
DE$eimpcnt[DE$eimpcnt==7]<-NA
DE$eimpcnt[DE$eimpcnt==8]<-NA
DE$eimpcnt[DE$eimpcnt==9]<-NA

#14 Richer

DE$eimrcnt
DE$eimrcnt[DE$eimrcnt== 7]<-NA
DE$eimrcnt[DE$eimrcnt== 8]<-NA
DE$eimrcnt[DE$eimrcnt== 9]<-NA


#15. poor

DE$imrsprc
DE$imrsprc[DE$imrsprc==7]<-NA
DE$imrsprc[DE$imrsprc==8]<-NA
DE$imrsprc[DE$imrsprc==9]<-NA

#16 Qualify: Education

DE$qfimedu
DE$qfimedu[DE$qfimedu ==77]<-NA
DE$qfimedu[DE$qfimedu ==88]<-NA
DE$qfimedu[DE$qfimedu ==99]<-NA


#17 Qaulify: Work Skill

DE$qfimwsk
DE$qfimwsk[DE$qfimwsk==77]<-NA
DE$qfimwsk[DE$qfimwsk==88]<-NA
DE$qfimwsk[DE$qfimwsk==99]<-NA


######################################################################


############################

  Run the original model

############################



## We ramdonly choose 500 samples

set.seed(123)
DE2 <- DE[sample(nrow(DE), 500), ]

original<-'

f1=~ ipeqopt + ipudrst + impenv + iphlppl + iplylfr

f2=~ ipmodst + imptrad + ipfrule + ipbhprp + impsafe + ipstrgv

f3=~imdfetn + eimpcnt + eimrcnt + imrsprc

f4=~ qfimedu + qfimwsk

#### missing f2~~f3, and f3~~f4

f1=~f4

f2=~f4
f1~~f2

f1=~f3


 '


## we need to include f1=~f3 to activate it

fit_original <- cfa(model = original, data=DE2, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")

fitMeasures(fit_original,c("chisq", "df", "pvalue", "nfi", "tli", "cfi", "rmsea"))

modindices(fit_original, minimum.value = 0, sort = TRUE)[1:10,]



newpar='

f1	~~	f3
f2	~~	f3
ipfrule	~~	ipbhprp
f3	~~	f4
f3	=~	ipudrst
ipeqopt	~~	imptrad
f2	=~	eimrcnt
f1	=~	ipfrule
impenv	~~	imptrad
f2	=~	iphlppl
 '

#### Univariate LM Test statitics

lavTestScore(fit_original, add=newpar)





#######################.  Bootstrap LM Test



newpar= ' 

f3~~f1
f3~~f2
ipfrule~~ipbhprp
f3~~f4
f3=~ipudrst
ipeqopt~~imptrad
f2=~eimrcnt
f1=~ipfrule
impenv~~imptrad
f2=~iphlppl'



set.seed(500)
Data<-DE2

bin.1000 <- matrix(NA, nrow = 500, ncol = 2)
for (i in 1:500) {
  
  boot.idx <- sample.int(nrow(Data), replace = TRUE)
  Data.boot <- Data[boot.idx,]
  
  fit2 <- sem(model = original, data = Data.boot, int.ov.free = FALSE, std.lv = FALSE, estimator = "ML")
  
  bin.1000[i, 1] <- lavTestScore(fit2, add = newpar)$uni[1, 4] # Don't forget to change from 1-10 manually
  bin.1000[i, 2] <- lavTestScore(fit2, add = newpar)$uni[1, 6] # Don't forget to change from 1-10 manually
}

t1 <- mean(bin.1000[, 1], na.rm = TRUE)
t2 <- mean(bin.1000[, 2], na.rm = TRUE)

t1
t2




##########################.   Bootstrap Wald Test


### in the following model, I didn't use f2~~f3, I insisted using f2=~f3



w_model<-'


f1=~ ipeqopt + ipudrst + impenv + iphlppl + iplylfr + b8*ipfrule

f2=~ ipmodst +  imptrad + ipfrule + ipbhprp + impsafe + + ipstrgv + b7*eimrcnt + b10*iphlppl

f3=~imdfetn + eimpcnt + eimrcnt + imrsprc + b5*ipudrst

f4=~ qfimedu + qfimwsk + b9*imrsprc

#### missing f1=~f3, and f3~~f4

f1=~f4
f1~~f2
f2=~f4



f1 ~~b1*f3
f2 =~b2*f3
ipfrule ~~b3*ipbhprp
f3 ~~ b4*f4
ipeqopt~~b6*imptrad
ipmodst ~~b8*imrsprc
impenv~~b9*imptrad
 '



## change b1 to b6 manually, one by one. 

con='

b1==0

 '


set.seed(500)
Data<-DE2

bin.1000 <- matrix(NA, nrow = 500, ncol = 2)
for (i in 1:500) {
  
  boot.idx <- sample.int(nrow(Data), replace = TRUE)
  Data.boot <- Data[boot.idx,]
  
  fit2 <- sem(model = w_model, data=Data.boot, meanstructure=FALSE, likelihood = "wishart", estimator = "ML")
  
  bin.1000[i, 1] <- lavTestWald(fit2, constraints = con)$stat[1]
  bin.1000[i, 2] <- lavTestWald(fit2, constraints = con)$p.value[1]
  
}

w1 <- mean(bin.1000[, 1], na.rm = TRUE)
w2 <- mean(bin.1000[, 2], na.rm = TRUE)

w1
w2




################.  The End

################. All results match the manuscript, 12/7/24




